library(spatstat)
## Loading required package: spatstat.data
## Loading required package: spatstat.geom
## spatstat.geom 3.2-9
## Loading required package: spatstat.random
## spatstat.random 3.2-3
## Loading required package: spatstat.explore
## Loading required package: nlme
## spatstat.explore 3.2-7
## Loading required package: spatstat.model
## Loading required package: rpart
## spatstat.model 3.2-11
## Loading required package: spatstat.linnet
## spatstat.linnet 3.1-5
##
## spatstat 3.0-8
## For an introduction to spatstat, type 'beginner'
Intensidad promedio estimada se puede calcular de dos formas diferentes
npoints(swedishpines)/area(Window(swedishpines))
## [1] 0.007395833
intensity.ppp(swedishpines)
## [1] 0.007395833
intensity(swedishpines)
## [1] 0.007395833
Más detalles sobre el proceso puntual se pueden obtener al realizar
summary(swedishpines)
## Planar point pattern: 71 points
## Average intensity 0.007395833 points per square unit (one unit = 0.1 metres)
##
## Coordinates are integers
## i.e. rounded to the nearest unit (one unit = 0.1 metres)
##
## Window: rectangle = [0, 96] x [0, 100] units
## Window area = 9600 square units
## Unit of length: 0.1 metres
Si existen marcas, se calculará la intensidad para cada tipo
plot(amacrine)
intensity.ppp(amacrine)
## off on
## 88.68302 94.92830
intensity(amacrine)
## off on
## 88.68302 94.92830
quadratcount(bei)
## x
## y [0,200) [200,400) [400,600) [600,800) [800,1e+03]
## [400,500] 271 401 100 32 215
## [300,400) 180 217 104 34 149
## [200,300) 334 4 26 33 172
## [100,200) 172 43 26 144 83
## [0,100) 146 89 234 338 57
bei |> quadratcount( nx = 10, ny=10) |> plot(main = "Quadratcount")
Si queremos calcular la intensidad de acuerdo al conteo obtenido por el quadratcount (puntos dividido para el área), se puede usar la función intensity sobre un cuadrante.
bei |> quadratcount( nx = 10, ny=10) |> intensity(, image=TRUE) |>
plot(main = "Intensidad respecto a quadratcount")
Suponiendo una distribución de Poisson, la hipótesis nula asume que el proceso es un CRS.
qt <- unmark(amacrine) |> quadrat.test(7,7)
qt
##
## Chi-squared test of CSR using quadrat counts
##
## data: unmark(amacrine)
## X2 = 9.6667, df = 48, p-value = 8.386e-10
## alternative hypothesis: two.sided
##
## Quadrats: 7 by 7 grid of tiles
plot(unmark(amacrine))
plot(qt, add = TRUE)
Comparar con la \(K\) función
env <- envelope(unmark(amacrine), Kest, nsim=19, rank=1, global=TRUE)
## Generating 19 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,
## 19.
##
## Done.
env
## Simultaneous critical envelopes for K(r)
## and observed value for 'unmark(amacrine)'
## Edge correction: "iso"
## Obtained from 19 simulations of CSR
## Envelope based on maximum deviation of K(r) from null value for CSR (known
## exactly)
## Alternative: two.sided
## Significance level of simultaneous Monte Carlo test: 1/20 = 0.05
## ...........................................................
## Math.label Description
## r r distance argument r
## obs hat(K)[obs](r) observed value of K(r) for data pattern
## theo K[theo](r) theoretical value of K(r) for CSR
## lo hat(K)[lo](r) lower critical boundary for K(r)
## hi hat(K)[hi](r) upper critical boundary for K(r)
## ...........................................................
## Default plot formula: .~r
## where "." stands for 'obs', 'theo', 'hi', 'lo'
## Columns 'lo' and 'hi' will be plotted as shading (by default)
## Recommended range of argument r: [0, 0.25]
## Available range of argument r: [0, 0.25]
## Unit of length: 662 microns
plot(env)
bei |>
envelope(Kest, nsim=19, rank=1, global=TRUE) |>
plot(main = "K para base bei")
## Generating 19 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,
## 19.
##
## Done.
Se asume que la función de intensidad \(\lambda(u)\) depende de la ubicación \(s\) en el espacio.
Los estimadores para la intensidad son:
\[ \tilde{\lambda}^{(0)}(u) = \sum_{i=1}^n \kappa(u-x_i)\]
Correlacionada uniformemente \[ \tilde{\lambda}^{(U)}(u) = \dfrac{1}{e(u)}\sum_{i=1}^n \kappa(u-x_i)\]
Corrección de Diggle \[ \tilde{\lambda}^{(U)}(u) = \sum_{i=1}^n \dfrac{1}{e(x_i)} \kappa(u-x_i)\]
para cualquier ubicación espacial \(u\) dentro de la ventana \(W\), donde \[ e(u) = \int_W \kappa(u-v) \]
y \(\kappa(u)\) es una función de densidad, es decir, \[ \kappa(u) \leq 0,\quad \forall u \in W \\ \int_{\mathbb{R}^2} \kappa(u) du = 1 \]
Observación
Una selección usual para el kernel es considerar una distribución normal.
La desviación estándar del kernel es el parámetro de suavizamiento o smoothing bandwidth. Un mayor valor de banda (bandwidth) incrementa el sesgo pero disminuye la varianza.
Teorema (Formula de Campbell)
Suponga \(f\) una función de ubicaciones espaciales \(u\), y considere la suma aleatoria \[ T = \sum_{i} f(x_i) \]
de los puntos \(x_i\) del proceso puntual \(X\). La formula de Cambell indica que \[ \mathbb{E}(T) = \mathbb{E}\left( \sum_{i} f(x_i) \right) = \int_{\mathbb{R}^2} f(u) \lambda(u) du \] donde \(\lambda\) es la función de intensidad de \(X\).
!!! pag 169, en la que se compara el uso de los estimadores.
Por defecto, se calcula el estimador uniforme del kernel, pero se puede indicar el estimador con la corrección de Diggle
amacrine |> density.ppp(diggle = TRUE) |>
plot(main = "Estimaición de la intensidad")
amacrine |> density.ppp(diggle = TRUE) |>
contour(main = "Estimaición de la intensidad")
amacrine |> density.ppp(diggle = TRUE) |>
persp(main = "Estimaición de la intensidad")
El valor de la banda indica el nivel de suavizamiento que tendrá la intensidad. Si no se especifica la banda, se toma por defecto un octavo de la longitud del lado más corto de la ventana considerada. Sin embargo, se puede seleccionar el valor optimo considerando por ejemplo
b <- bw.diggle(amacrine)
## Warning: Berman-Diggle Cross-Validation criterion was minimised at right-hand
## end of interval [0, 0.0624]; use argument 'hmax' to specify a wider interval
## for bandwidth 'sigma'
b
## sigma
## 0.06237769
plot(b)
amacrine |> density.ppp(diggle = TRUE, sigma = b) |>
plot(main = "Estimaición de la intensidad con valor óptimo")
La banda puede variar bastante dependiendo del método que se utilice, donde cada uno tiene supuesto que pueden no adaptarse al modelo estudiado.
bw.ppl asume un proceso de Poisson
inhomogéneo.
bw.diggle asume un proceso de Cox, que forma más
grupos o clustes clustered que un proceso de Poisson.
bw.scott
En todos los casos, hemos supuesto que tanto como la función de
suavizamiento como la banda es la misma para todo el proceso. Eso puede
mejorarse si se consideran métodos adaptativos. Referencia el paquete
sparr en R.
Aquà f representa la fraacción de números del patrón
puntual tomados de manera aleatoria y usados para realizar la partición
de Dirichlet.
amacrine |> adaptive.density(f=1) |> plot()
amacrine |> adaptive.density(f=0.1, nrep = 30) |> plot()
## Computing 30 intensity estimates...1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
## 21, 22, 23, 24, 25, 26, 27, 28, 29,
## 30.
## Done.
Otra opción es utilizar el vecino más cercano
amacrine |> nndensity(k = 10) |> plot()
Intensidad homogénea: \(\lambda(u) = \lambda\)
Intensidad homogénea por regiones \(\lambda(u) = b_j\) para todo \(u \in B_j\).
Intensidad proporcional a un valor base: \(\lambda(u) = \theta b(u)\)
Modelo exponencial: \(\lambda(u) = \kappa e^{\beta Z(u)} = \exp(\alpha + \beta(Z(u))\) donde \(Z(u)\) es una covariable espacial y \(\kappa\) y \(\beta\) son parámetros.
Modelo de incidencia exponencial: \(\lambda(u) = b(u) \exp(\alpha + \beta Z(u))\)
Modelo log lineal: \(\lambda_{\theta}(u) = \exp(B(u) + \boldsymbol{\theta}^\intercal \boldsymbol{Z}(u) ) = \exp(B(u) + \theta_1 Z_1(u) + \theta_2 Z_2(u) + \ldots + \theta_p Z_p(u) )\)
Ubicación de árboles en una zona. Para el modelo log lineal
\[\lambda(u) = \exp(\beta_0 + \beta_1 S(u))\] donde \(S(u)\) es la pendiente en la ubicación \(u\).
fit <- ppm(bei ~ grad, data=bei.extra)
fit
## Nonstationary Poisson process
## Fitted to point pattern dataset 'bei'
##
## Log intensity: ~grad
##
## Fitted trend coefficients:
## (Intercept) grad
## -5.391053 5.026710
##
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## (Intercept) -5.391053 0.03001787 -5.449887 -5.332219 *** -179.5948
## grad 5.026710 0.24534296 4.545847 5.507573 *** 20.4885
de donde \(\beta_0 \approx -5.39\) y \(\beta_1 \approx 5.03\).
Podemos observar el efecto que tiene la covariable en la intensidad
plot(effectfun(fit, "grad", se.fit=TRUE))
Otro modelo a considerar serÃa
\[\lambda(u) = \exp(\mu + \alpha_{E(u)} + \alpha_{G(u)} + \gamma_{E(u),G(u)})\] es decir,
fit <- ppm(bei ~ elev * grad, data=bei.extra)
fit
## Nonstationary Poisson process
## Fitted to point pattern dataset 'bei'
##
## Log intensity: ~elev * grad
##
## Fitted trend coefficients:
## (Intercept) elev grad elev:grad
## -4.403426761 -0.007024363 -36.500458701 0.292813497
##
## Estimate S.E. CI95.lo CI95.hi Ztest
## (Intercept) -4.403426761 0.614763210 -5.60834051 -3.198513010 ***
## elev -0.007024363 0.004192219 -0.01524096 0.001192235
## grad -36.500458701 5.261456400 -46.81272375 -26.188193651 ***
## elev:grad 0.292813497 0.036257125 0.22175084 0.363876155 ***
## Zval
## (Intercept) -7.162801
## elev -1.675572
## grad -6.937330
## elev:grad 8.076026
Podemos graficar también la estimación obtenida
plot(fit, how="image", se=FALSE)
Confidence intervals for the parameters
confint(fit, level=0.95)
## 2.5 % 97.5 %
## (Intercept) -5.60834051 -3.198513010
## elev -0.01524096 0.001192235
## grad -46.81272375 -26.188193651
## elev:grad 0.22175084 0.363876155
Depositos de oro en Murchison, Australia. Primero pasamos los datos a kilómetros y calculamos la distancia a las fallas para cada punto.
mur <- solapply(murchison, rescale, s=1000, unitname="km")
dfault <- with(mur,distfun(faults))
estimados el modelo
\[\lambda(u) = \exp(\beta_0 + \beta_1 d(u))\] donde \(d(u)\) es la distancia más corta desde \(u\) a la falla geológica más cercana.
fit <- ppm(gold ~ dfault, data=mur)
fit
## Nonstationary Poisson process
## Fitted to point pattern dataset 'gold'
##
## Log intensity: ~dfault
##
## Fitted trend coefficients:
## (Intercept) dfault
## -4.3412775 -0.2607664
##
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## (Intercept) -4.3412775 0.08556260 -4.5089771 -4.1735779 *** -50.73802
## dfault -0.2607664 0.02018789 -0.3003339 -0.2211988 *** -12.91697
de donde \(\beta_0 \approx -4.34\) y \(\beta_1 \approx -0.26\).
plot(effectfun(fit, "dfault", se.fit=TRUE))
Se compara con un modelos CRS y en este caso se rechaza.
fit0 <- ppm(bei ~ 1)
fit1 <- ppm(bei ~ grad, data=bei.extra)
anova(fit0, fit1, test="LR")
## Analysis of Deviance Table
##
## Model 1: ~1 Poisson
## Model 2: ~grad Poisson
## Npar Df Deviance Pr(>Chi)
## 1 1
## 2 2 1 383.11 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Se estudia la inclusión o no de ciertas variables
fit2e <- ppm(bei ~ polynom(elev, 2), data=bei.extra)
fit2e1g <- update(fit2e, . ~ . + grad)
anova(fit2e, fit2e1g, test="LR")
## Analysis of Deviance Table
##
## Model 1: ~elev + I(elev^2) Poisson
## Model 2: ~elev + I(elev^2) + grad Poisson
## Npar Df Deviance Pr(>Chi)
## 1 3
## 2 4 1 419.24 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Preferimos el modelo con el AIC más pequeño.
fitprop <- ppm(bei ~ offset(log(grad)),data=bei.extra)
fitnull <- ppm(bei ~1)
AIC(fitprop)
## [1] 42494.42
AIC(fitnull)
## [1] 42763.92
diagnose.ppm(fit2e)
## Model diagnostics (raw residuals)
## Diagnostics available:
## four-panel plot
## mark plot
## smoothed residual field
## x cumulative residuals
## y cumulative residuals
## sum of all residuals
## sum of raw residuals in entire window = 1.62e-09
## area of entire window = 5e+05
## quadrature area = 5e+05
## range of smoothed field = [-0.005641, 0.006091]
fit <- ppm(bei ~ polynom(grad, elev, 2), data=bei.extra)
lamhat <- predict(fit)
lamB <- predict(fit, locations=bei)
predict(fit) |> plot()
X <- simulate(fitprop);plot(X[[1]])
# Pesos
vols <- with(marks(finpines),
(pi/12) * height * (diameter/100)^2)
Dvol <- density(finpines, weights=vols, sigma=bw.ppl)
plot(Dvol)